// M Orsi - May 2007. This program computes the force autocorrelation function needed to predict the diffusion
// coefficient through the x-constraint algorithm.
// Input data: a file called "constraintForce.dat" containing in each line the instantaneous value of the constraint force
// Input parameters: file "tAcf.in" containing the desired length, in steps, of the facf
// Output: autocorrelation function of the random force
// Main ref: Allen and Tildesley, pp. 186-187

// Compilation: gcc -Wall -O -o facf facf.c -lm

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define	kB_J_K		1.3806505e-23 /* [J / K] */
#define AVOGAD		6.0221367e23 /* [mol^(-1)] */
#define kB_J_K_mol      (kB_J_K * AVOGAD) // [ J / (K mol) ]
#define kB_kJ_K_mol      (kB_J_K_mol / 1000) // [ J / (K mol) ]
#define T                303  // temperature [K]
#define kT           ( kB_kJ_K_mol * T ) /* [ kJ / mol) ] */
#define kT2           ( kT * kT ) 

#define AllocMem(a, n, t)  a = (t *) malloc ((n) * sizeof (t))
#define Min( x1, x2 )  ( ( ( x1 ) < ( x2 ) ) ? ( x1 ) : ( x2 ) )

// global variables
double *f, // force array
  *facf,   // force autocorralation function (FACF) array
  *xi,   // friction coefficient (FACF) array
  //  *norm, // auxiliary array
  intXi,  // time-integral of xi, ie., static friction coeff
  intFacf; // time-integral of FACF
int tRun;  // # force data = #timeSteps of the run 
int tAcf;  // facf length

const double timeStep_ns = 20e-6;      // [ns] this value is actually 20 fs, but it's better to convert it to ns
const double timeStep_ps = 20e-3;      // [ns] this value is actually 20 fs, but it's better to convert it to ns
/* const double acfLengthFactor = .002; // this defines the acf time-window in relation to the total data-collection time */
/*                                    // Large values will yield long correlation functions but poor statistics, and vice-versa  */

void AllocArrays()
{
  FILE *fPtr;
  char buff[80];
  // count data
  if ( ( fPtr = fopen( "constraintForce.dat", "r" ) ) == NULL)
    printf( "File constraintForce.dat could not be opened\n" );
  else {
    tRun = -1;
    while ( !feof( fPtr ) ) {
      fgets( buff, 80, fPtr );
      //printf("tRun = %d\n", tRun);
      tRun++;
    }  
    fclose( fPtr );
  }
  if ( ( fPtr = fopen( "tAcf.in", "r" ) ) == NULL) printf( "File tAcf.in could not be opened\n" );
  else {
    fscanf( fPtr, "%d", &tAcf );
    fclose( fPtr );
  }
  AllocMem( f, tRun, double );
  AllocMem( facf, tAcf + 1, double ); // +1 is needed, see p.186 AT third line of code
  AllocMem( xi, tAcf + 1, double ); // +1 is needed, see p.186 AT third line of code
}

void  LoadForceData()
{
  FILE *fPtr;
  double force;
  int i = 0;
  
  if ( ( fPtr = fopen( "constraintForce.dat", "r" ) ) == NULL)
    printf( "File constraintForce.dat could not be opened\n" );
  else {
    while ( !feof( fPtr ) ) {
      fscanf( fPtr, "%lf", &force );
      f[ i ] = force;
      i++;
    }
  }
}

void  ComputeRandomForce()
{
  double avgForce = 0.; // time-averaged Fcz
  int i;
  
  for ( i = 0; i < tRun; i++ ) avgForce +=  f[i];
  avgForce /= tRun;
  printf( "% 6.3f kJ / (mol nm) <- Average constraint force Fzc\n", avgForce );
  
  // Calculate random forces, i.e., deltaF = instantaneousF - avgF
  for ( i = 0; i < tRun; i++ ) f[i] -= avgForce;

  // printing for debugging purpose
  FILE *fPtr;
  if ( ( fPtr = fopen( "rndmForce.dat", "w" ) ) == NULL )
    printf( "File rndmForce.dat could not be opened\n" );
  else { 
    for ( i = 0; i <= tRun; i++ ) fprintf( fPtr, "%f\n", f[ i ] );
    fclose( fPtr );
  }
}

void  ComputeFacf() // ref: allen-tildesley p.187
{
  int t, tMax, t0;
  double ftf0;

  //  printf( "Computing FACF... " );
  for ( t = 0; t <= tAcf; t++ ) {
    tMax = tRun - t;
    ftf0 = 0.;
    for ( t0 = 1; t0 <= tMax; t0++ ) ftf0 += f[ t0 ] * f[ t0 + t ];
    facf[ t ] = ftf0 / tMax;
  }
}

void  PrintFacf()
{
  int t;
  FILE *fPtr;
  char facfFileName[ 12 ];
  char time_ps[ 6 ]; // string containing timeNow, e.g. "100", in ps

  sprintf( time_ps, "%5.5d", ( int ) ( 0.5 + tAcf * timeStep_ps ) );
  facfFileName[ 0 ] = '\0';
  strcat( facfFileName, time_ps );
  strcat( facfFileName, "ps.facf" );

  if ( ( fPtr = fopen( facfFileName, "w" ) ) == NULL ) printf( "File .facf could not be opened\n" );
  else { 
    for ( t = 0; t <= tAcf; t++ ) { 
      fprintf( fPtr, "%f\t%e", t * timeStep_ps, facf[ t ] * 2.76e-24 ); // [ ps ]  [ N^2 ]
      //      fprintf( fPtr, "\t%f", facf[ t ] ); //                               [ ps ]  [ kJ^2 / ( mol^2 nm^2 )] 
      fprintf( fPtr, "\n" );
    }
    fclose( fPtr );
  }
}

void IntegrateFacf() // ref: Rapaport p. 491
{
  int i;
  
  intFacf = .5 * ( facf[ 0 ] + facf[ tAcf ] ); 
  for ( i = 1; i < tAcf; i++ ) intFacf += facf[ i ];
  intFacf *= timeStep_ns; // time-scaling
  //  printf( "Time-integral of FACF: %6.3f  kJ^2 ns / ( mol^2 nm^2 )\n", intFacf );
}

void ComputeFrictionCoefficient()
{
  int i;
  FILE *fPtr;

  char xiFileName[ 12 ];
  char time_ps[ 6 ]; // string containing timeNow, e.g. "100", in ps

  sprintf( time_ps, "%5.5d", ( int ) ( 0.5 + tAcf * timeStep_ps ) );
  xiFileName[ 0 ] = '\0';
  strcat( xiFileName, time_ps );
  strcat( xiFileName, "ps.xi" );

  if ( ( fPtr = fopen( xiFileName, "w" ) ) == NULL ) printf( "File .xi could not be opened\n" );
  else { 
    for ( i = 0; i <= tAcf; i++ ) {
      xi[ i ] = facf[ i ] / kT ;
      fprintf( fPtr, "%f\t%f\n", i*timeStep_ps, xi[ i ] ); // [ ps ]   [ kJ / ( mol nm^2 )]
    }

    // compute integral of the friction coefficient, i.e., static friction coefficient
    intXi = .5 * ( xi[ 0 ] + xi[ tAcf ] );
    for ( i = 1; i < tAcf; i++ ) intXi += xi[ i ];
    intXi *= timeStep_ns; // time-scaling [ns]
    //    printf( "Static friction coefficient xi = %6.3f kJ * ns / (mol nm^2)\n", intXi ); // [ kJ ns / (mol nm^2) ]
    fclose( fPtr );
  }
}

void ComputeDiffusionCoefficient()
{
  double Dz; // diffusion coefficient
  
  //Dz = kT2 / intFacf; // units!???? [ nm^2 / ns ] = [ m^2 / s ]
  //printf ("kT = %6.3g kJ / mol\n", kT);
  Dz = kT / intXi; // [ nm^2 / ns ]
  //  printf( "Diffusion coefficient = % 5.3E nm^2 / ns\n", Dz );
  printf( "% 5.3E cm^2 / s <- Diffusion coefficient\n", Dz * 1e-5 );
  //  printf( "Diffusion coefficient = % 5.3E nm^2 / microsec\n", Dz * 1e3 );
}

int main()
{
  AllocArrays();
  LoadForceData();
  ComputeRandomForce();
  ComputeFacf();
  PrintFacf();
  IntegrateFacf();
  ComputeFrictionCoefficient();
  ComputeDiffusionCoefficient();
  printf( "tRun = %d = # force data = # lines in the constraintForce.dat file\n", tRun);
  printf( "tAcf = %.1f ps\n", (tAcf + 1) * timeStep_ps );
  return 0;
}

/*   // Alternative facf calculation from allen-tildesley p.186 
    This one should be faster, but it isn't! So keep the above code which is even clearer */
/*   for ( t = 0; t <= tAcf; t++ ) facf[ t ] = 0; */
/*   for ( t = 0; t <= tAcf; t++ ) norm[ t ] = 0; */

/*   double f0; */
/*   int tt0Max, tt0; */

/*   for ( t0 = 1; t0 <= tRun; t0++ ) { */
/*     f0 = f[ t0 ]; */
/*     tt0Max = Min( tRun, t0 + tAcf ); */
/*     for ( tt0 = t0; tt0 <= tt0Max; tt0++ ) { */
/*       t = tt0 - t0; */
/*       facf[ t ] += f0 * f[ tt0 ]; */
/*       norm[ t ] += 1.; */
/*     } */
/*   } */
/*   for ( t = 0; t <= tAcf; t++ ) facf[ t ] /= norm[ t ]; */
